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ABSTRACT 

A new technique for constructing "computational molecules" 
for linear finite difference operators is developed. The basic 
approach is one of approximating a two dimensional surface with 
a geometrically consistent interpolating polynomial of degree 
four or five. The desired finite difference operator is then 
developed from the polynomial. The resulting molecules are 
geometrically consistent and may be used to solve boundary 
value problems without the use of fictitious points. 

Molecules for the biharmonic operator with various 
boundaxy conditions included are presented in this paper, as 
well as molecules representing the boundary conditions for 
shear and moment along the free edge of a plate. 

The integrity of the molecules presented is proven by 
comparison of solutions for flat plate bending problems by 
finite difference with exact solutions from the literature. 


Convergence plots for each problem are also presented. 
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NOMENCLATURE 


Differentiation of some function A with respect 
to variables (a) and (b). 


Constants for an interpolating polynomial of degree 
if Owls 


heats} 3 
Flexural rigidity: — ET 
iene) 
Constants for an interpolating polynomial of degree. 
Young's modulus of elasticity. 
Error of approximation. 


Expression for variable terms of an interpolating 
polynomial of degree four. 


Expression for variable terms of an interpolating 
polynomial of degree five. 

Mesh size spacing. 

Length or width. 


Memenesln tne normal direction along the edge of a 
flat plate. 


Number of divisions on one side of a square finite 
difference mesh. 


Normal Co-ordinate. 

Applied pressure 

Interpolating polynomial. 

Values of polynomial at various nodal points. 
Tangential co-ordinate. 

Transverse shear along a free edge of a plate. 
Horizontal co-ordinate. 

Vertical co-ordinate 

Plate deflection. 


Pousson s rat1o 





MNaeeieemoratton Lor column vector. 
Matai netatlon tor row vector. 


Square or rectangular matrix. 


Harmonic operator. 


Biharmonic operator. 
Thickness of a thin plate. 
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I. INTRODUCTION 

A. GENERAL DESCRIPTION 

Solving partial differential equations and boundary value 
problems by approximate numerical methods is not new. Such 
numerical solutions have become very popular in the past twenty 
years due to the development of the inate speed digital computer. 

One widely used numerical method for solving partial 
differential equations is that of finite differences. The 
solution of equations by this method leads to the development 
of “computational molecules". Various techniques have been 
used to develope such molecules and numerous molecules are 
available in the literature (Refs. 1,2,3,&4). 

Often these molecules tend to be geometrically inconsistent. 
This inconsistency results from combining molecules for the 
elementary differential operators, developed from different 
interpolating polynomials, in order to obtain a more complex 
operator. For instance, the molecule for the biharmonic operator 
as given by Salvadori and Baron (Ref. 1s} utilizes two different 
interpolating formulas. The fourth derivatives are obtained 
from a polynomial of degree four. The mixed partial derivative 
1s obtained from a bi-quadratic interpolating surface. 

Furthermore, the application of these molecules in the 
solution of boundary value problems requires the use of ficti- 
tious nodal points outside the physical boundaries of the 
problem. Relationships between these fictitious nodal points 
and those nodal points representing the physical problem must 
be pre-determined. When the fictitious points cannot be 


determined uniquely, the method fails. 
7 





Semetoeming the above disadvantages of the classical finite 
difference operators found in the literature, the objective 
of the author's research was to develope a new technique for 
constructing linear finite difference operators. This new 
technigue utilizes one consistent form of an interpolating 
polynomial to approximate a two dimensional surface from which 
molecules for all operators are constructed. Furthermore, 
basic operators incorporating physical boundary conditions 
are constructed using this technique. The use of such operators 
alleviates the need for fictitious nodal points. 

The technique proposed in this research is applicable 
for the construction of any finite difference operator. How- 
ever, the field of investigation for the purpose of this 
research was limited to the construction of those operators 
required for the solution of plate bending problems with various 
boundary conditions. The solutions of such problems will be 
used to justify the integrity of the constructed finite dif- 


ference operators. 


re BASIC TECHNIQUE 

The basic technique employed to construct the molecules 
for a prescribed linear finite difference operator is one of 
approximating a two dimensional surface by means of a biquad- 
PaAetewor O1quIntte interpolating polynomial. The constants 
for such an interpolating polynomial are determined for each 
desired surface, which may include specified boundary condi- 


Dione wmeuSincgwaltinite @irference grid. The polynomial is 





then operated on with the desired partial differential 
operator to be approximated. 

The final result is a numerical approximation of the 
Operator in terms of a computational molecule, grid size, 
nodal deflections, and nodal boundary condition values. 

All calculations required for the construction of these 
molecules were performed in double precision using an 


Pey—-35600 computer. 
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fio cove ON OF DIFFERENTIAL OPERATORS 


A. INTERPOLATING POLYNOMIALS 
The interpolating polynomials used to approximate a two 
dimensional surface for the development of a finite difference 


biharmonic operator in this work were either 


P(x,y) = Cyy?x*+ Coy?x "+ Cay?x 44 Cyyx*+ Cray 


+ ae eee aso ome Leas 


y4x24C. y3x24+C_, y2x24C_ ,yx24+C__x2 


Saiki 2 eS 14 al) 


4 3 2 
+C1 ¢6Y Xx +C)9Y xX +C1 gy 6 +C) gYX +C5 9X 


+Co1y4 +Co9y? +C53y7 +C94Y ee 1G (2.1) 
and in matrix notation 
(1) 
eB  (xry)? {€3} 5-1 2,...,25 (2. 1a) 
or 
P(x,y) = Daa Doy*x?+ Degas Dysy?x?+ Deyx>+ Dax 


55-4 4,4 3,4 25.4 4 4 
+ Doy xo Dey x + Doy Xx +D, 9Y°x +D,7Y* +Dy 5% 
Dy 3 453 3y3 23 3 x3 
+D, 3Y~x +D, AY x +D)5Y~X +D1 ¢Y x +D)7Yx +D) 9x 
Dene 4,2 oe Deve Z 2 
+D) oY x +Do oY Xx +D51Y Xx +D55Y Xx +D5 3Y%X +Do 4X 
oD. 3s +D yox +D y 3x +D yx tay Xn DX 
2 26 2/7 28 Zo 30 


Cees FD4y- +D 


ane 32Y Bay ae 
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and in matrix notation 
P(x,y) = (G6) (x,y) ) (Dee 
fine eee 3 6 (2.2a) 

Polynomials (2.1) (2.2) will remain polynomials of degree 
four or five respectively for both constant values of x and 
y. All partial derivatives with respect to x and y of order 
four exist. Either polynomial will consistently approximate 
a two dimensional surface. 

In certain cases, a fifth order polynomial such as (2.2) 
which will approximate a thirty-six nodal point grid surface 
is required in the construction of the biharmonic operator 
molecule with included boundary conditions. This requirement 
results from the need for keeping the error of approximation 
for a given aatecuie @£ an order equal to or greater than the 
square of the grid size spacing. The subject of approximation 


Error will be further discussed later. 


b- BASIC FINITE DIFFERENCE BIHARMONIC OPERATOR 

For the solution of a problem governed by a partial 
differential equation, in which the biharmonic operator appears, 
a numerical approximation of the biharmonic operator without 
boundary conditions must be constructed first. To construct 
the molecule for such an operator ~olynomial (2.1) was used to 
approximate a two dimensional s se defined by the ordinates 
at the nodes of a grid shown in fase 11.) | 

To determine the unknown constants for the interpolating 


polynomial, the polynomial was evaluated at each of the 
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twenty-five nodal points in terms of the co-ordinate values 
of x and y. In matrix notation, the resulting system of 


equations are expressed as 


eee) = ed a eT G23) 
where 
es {P(xisYi)Fsiy 201, 25 


[M} ( F]) (x3,¥))) 


( EF) (x5,y9)) 


Gace.) 
ee a | 


Solving equation (2.3) for the unknown coefficients, yields 


: = ml ; 
ici} [IM] {Pisa a 6.25 (2.4) 


with which the interpolating polynomial is rewritten as 


eee 


(i) 
Creeley) [M] ey Oo 25a 2 5) 


Pixy) =" 


Operating on equation (2.5) with the biharmonic operator 


results in the following equation 


alee 


4 eyes a) =1 
ee ee eee) IM) 1Piy g 2s (2.6) 


The quantity ge FF) (x y)) is now evaluated for a particular 


nodal point say m(x=2h,y=2h). The resulting quantity 
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( ve F) (2h, 2n)) [mM] 71 after numerical evaluation yields a 
five by five computational molecule for the finite difference 
biharmonic operator at nodal point (m). 

The above molecule allows for the numerical approximation 
of the biharmonic operator in terms of known constants and 
nodal deflections. Evaluation of a computational molecule 
for a nodal point other than the center one is performed in 
a Similar manner except that polynomial (2.2) and a six by 
Six grid is used in order to achieve an acceptable bound on 
the error of approximation. 

Cc. FINITE DIFFERENCE BIHARMONIC OPERATORS WITH BOUNDARY 

CenOTTTONS SPECIFIED 

To complete the solution of a boundary value problem, the 
applicable boundary conditions must be considered. If geo- 
metric considerations alone are required for the boundary 
conditions, a general finite difference molecule for the 
biharmonic operator incorporating the desired boundary con- 
ditions may be developed. The technique uses polynomial (2.2) 
to approximate a surface having the desired boundary conditions. 


Vs Operators With First Partial Derivative Boundary 
Conditions 


HOGmenenCOncenuctLon Of a finite difference bihar- 
monic operator with a first partial with respect to y boundary 
condition, polynomial (2.2) is used to approximate a surface 
such as the one represented by Figure (2) where the nodal 


points along the top edge are repeated to accommodate both 
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values of deflection and the specified boundary conditions 
mlong that edge. 

The unknown coefficients for the interpolating 
polynomial are determined by evaluating the polynomial at each 
of the thirty distinct nodal points and the first partial with 
respect to y of the polynomial at each of the nodal points 


along the top edge. The resulting equations are 


P(xy) = (6%) (x,y) ) w= 26 (2.7) 


P iy (XY) (6) (x,y) ) esi 36 (2.8) 


The complete system of equations becomes 


Tes = [My J {Dito 2,...,36 (2232) 
where 
Jb 
{P.} = (P(xpr¥5) | Phy Or ¥K) ? 
= cy) 
[My } 47” 5-1 2;...,30 


and solving for the unknown coefficients 


Vale aac 0 (2.10) 


> 





from which the interpolating polynomial may be written as 
P(x,y) = (G@)(x,y))  [M,17* {P,) 
eee sce 57 2 0 (2341) 
Operating on equation (2.11) with the biharmonic operator 
results in the following equation 


vipix,y) = (V4GG) (x,y)) tm] {23} 
eee sp SO 2 ee) 


The quantity a ee 1s evaluated for any 
desired nodal point within the grid and the resulting numerical 
evaluation of (y4e6™) (x,y) ) eal gives a Six by six molecule 
for the biharmonic operator with the prescribed boundary 
conditions. 

For first partial derivatives with respect to x and 
to y along two perpendicular edges, the respective biharmonic 
molecule is developed using polynomial (2.2) to approximate 
the surface defined by ordinates at the nodes of the grid 
given in Figure (3). The nodal points along two perpendicular 
edges are repeated to account for both the values of deflections 
and the boundary conditions at these nodal points. 


The coefficients for the interpolating polynomial are 


determined from the following relationships 


( G2) (x,y) ) {D. } 


P (x,y) 
1 i=1,2,...,36 (2s) 


( of) (x,y) ) fen 


Py (Xry) a 


MM 2 ts eeiged O (2.14) 
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Le 





II 


( . con {D; } 


P (x,y) 
f epee. 36 (2.15) 


(GO) Gy)? {D4} 


P CY) 
,XY Oe ae, o0 (2.16) 


The mixed second partial was used in the development in order 
to have sufficient conditions to evaluate the thirty six 
coefficients of the interpolating polynomial and to provide 
for continuity at the common corner nodal point. 


The complete system of equations becomes 


ee ee eee 
fas 6 7 30 (22) 
where 
7 q 
eer =n PCRS YG VIP y (RR rYx) IP Xp eV) IP ey OY) 


(sae 
Go" (a5 ry5)? Fe te 5 
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Solving for the constants, {D,;}, and operating on 
the resulting polynomial with the biharmonic operator yields 
vep(x,y) = ( V4G) (x,y) ) tMgi7t (P,)- 

Alegre eats Od (2 .ako) 
where (vig |) (x,y) ) [M,]~1 is evaluated for any desired 
nodal point. The result is a molecule for the biharmonic 
operator with the desired specified boundary conditions. 


2 Crcmceoiseitemsoccond Partial Derivats; '¢ Boundary 
Conditions 


A biharmonic operator for a region with the harmonic 
operator specified along a part of the boundary is often 
reguired for the solution of a boundary value problem. The 
required molecule for the biharmonic operator with such speci- 
fied boundary conditions is developed using polynomial (2.2) 
and the grid of Figure (2). The applicable equations for 


determining the constants of the polynomial are 


(64) (x,y)) {p:} 


ae 


Pix, y) 
Wiis ec «70 (2.19) 


Il 


v-p (x,y) = (Vv26 4) (x,y) ) {D2} 


pees eg oO 2210) 


from which the system of equations becomes 


{p.} = [M]} {D.} 
= : wo. 36 (2.1) 
where 
| 2 = 
Lae) = CP (xs ,y5)!V Ey) 


io 





i= 2; ><a Oe 


Solving for the unknown constants and operating on the 
resulting polynomial with the biharmonic operator yields the 
desired relationship for determining the biharmonic molecule. 
The relationship for the operator with the specified boundary 
eonaition is 
V{P (x,y) = (y4g%) (x,y)? [M3] ~* {D, } 
teres SO. e le ace) 

For boundary conditions requiring the second partial 
with respect to y along one edge and the second partial with 
respect to x along a perpendicular edge, polynomial (2.2) and 
the grid of Figure (3) is used in the development of the 
applicable biharmonic molecule. The equations required to 
determine the appropriate constants for the interpolating 
polynomial are 


(6) (x,y)) {Dz} 
ieee ee 7 SO (2 223)) 


Px, y) 


P (oy) 
yy 


f 


Ger iGay) 1D, | 
pey, Pee 7 oie OO (2.24) 
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ee 
Pex (ry) = (GS) Gy)? {Dy} 


aap? pe 36 (2.25) 

VP (x,y) = (yc!) (x,y)) {D; } 
=e Pet =. 36 (2226) 

P viery) = (6%?) (x,y)) {0} 
yen rXXY 2 eae 36 (Qua 
1YYX ails i=99 2) s,, ,26 (2220) 


Conditions expressed by equations (2.26), (2.27), and (2.28) 

are required to obtain sufficient independent equations with 
which to determine the thirty-six unknown constants. The 

values for these additional nodal boundary conditions are 
handled as additional unknowns. The resulting set of equations 
for determining the unknown constants of the required polynomial 


is 


tee 


il 
= 
> 
Pan 
oO 
- 
Lynd 
je 
il 


ei , 30 (2.29) 
where 


i | | 
, Jk 
zee 2G 


Re ee 35) 2 yn 36") 36) ” 
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( Gy 361 36)) 


tel Zar 
The final form of the desired finite difference 
biharmonic operator is 
vip(x,y) = (v46%) (x,y) [M)7) {D,) 
ee peereetig On here 3 0) 
Bi Operators With Mixed Boundary Conditions 
The development of a biharmonic molecule with boundary 


ConcmetOns Or the first partial with respect to y along one edge 


(a)tt should be noted that points 26, 35, and 36 all have 
the same geometrical co-ordinates. 
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Peete harmonic Condition along a perpendicular edge utilizes 
Meeerpolating polynomial (2.2) and the grid of Figure (3). The 
applicable equations for determining the required constants 


for the interpolating polynomial are 


P (x,y) = (G6) (x,y)) {D,} 
i=1,2,...,36 (remy) 

P (x,y) = (G6) (x,y))  {p.} 
'Y ry Aeon 36 (D3 on 

V*Pix,y) = (026%) (x,y)) {D,) 
j=1,2,...,36 (83) 
(1) 

P (x,y) = ¢(G Gay) {D.} 

gy XY le eee G (DAREN 


The resulting system of equations is 


{P.} = [M,] {D;} 
z en) en, 36 (2.35) 
where 
| | > ( 7 
Weed CP (xs IP Y Oyry,) IV P(Xor¥g) IP yo (%36 036)? 


Zs 





KaziGy 2 pes 7 OU 


—_ = fF 2 0° = fF  F— lS Ol Sl lll ESESE-lClelllClCellC 


Ty ise, we 
tice taakerormneOor Ene desired finite difference 


biharmonic operator is 


vip(x,y) = (v46™) (x,y))  [Mg]72 {Dy}. 


Hi | poe O01 eZee) 
4. Other Finite Difference Operators 

Certain physical boundary conditions include values 
of material properties in their evaluations. For example, in 
plate bending two such conditions are the equations for the 
moment and the shear along the free edge of a flat plate under 
some loading. The technique of combining the biharmonic 
operator and such boundary conditions into a single finite 
difference molecule is not practical in such cases. The reason 
being that in order to invert the matrices encountered in the 


development of the desired molecule, the material property 
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would have to be specified before the inversion of matrix [M]. 
This would limit the application of the molecule to a specific 
material. 

To prevent this limitation, a molecule for the 
desired boundary condition alone is developed using the con- 
Sistent interpolating polynomial (2.1). The technique is one 
of developing individual molecules for the partial derivatives 
involved in the boundary condition. The individual molecules 
are then combined with the applicable physical properties to 
form a general molecule which will represent the boundary 
@ondition. 

An example of this approach is illustrated in the 
development of molecules for the moment and the shear along 


the free edge of a flat plate. The governing equations are 


(Ref. 5) 

Mn = D(P nn (X,Y) to P etlx,y)) = Q (2.37) 
and 

ene M2) P ae. (x;y) = 0 (2.38) 


where n is the normal co-ordinate, t the tangential co-ordinate, 
mies) Poisson's ratio and D the flexural rigidity of the plate. 
Molecules for each of the derivatives in equations 
(2.37) and (2.38) are developed by using the interpolating 
polynomial in the form given by (2.5). Operating on the 
interpolating polynomial with the required differential oper- 


ators yields the following r«:lationships 
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Pp oy) = ( Ff) (x, y)) [M]~1 {p.} 
oe ae es 29) 


P (x,y) = €F@)(x,y)) [72 {P.} 
pel tt aeRO 5 We Aa) 
= 
P (x,y) = (FF) (x,y)) [m7 * {Pe} 
a ynnn ety ee Oma) 
P ttn (XY) = ( PAP) Gy)? (mM) 71 ee 


i=1,2,..-25 (2.42) 


Evaluation of the above relationships at a nodal point along 
the boundary of the grid given in Figure (1) yields the 
required molecule for each derivative. 

The molecules are then combined to form the following 


finite difference approximations for equations (2.37) and 


(230). 
= lt -1 
My = D(C PG) (x,y) EM to (PGP Gy) EN) {Py} 
C2 44) 
= L) ol. E (i) -1 
vy = o(C FO) (x,y)) im + (2-n) CPO) yy) Ea? (Py) 


(2.44) 


D. ERROR OF APPROXIMATION 

The molecules developed in this work are based on approxi- 
mating a two dimensional surface wit’ a polynomial. Before the 
results can be effectively used, a bound on the error entailed 


in this approximation is essential. 


26 





The error of approximation for a finite difference 


biharmonic operator for some nodal point (m) is defined as 


4 4 4 
0°P 0 °F ol 4 
: am a? ay? 'Y) (APPROX.) ~—- (2.45) 


The approximate value of the biharmonic operator at nodal point 
(m) is determined from the appropriate computational molecule 
and the nodal values associated with the molecule. 

An evaluation of the required nodal values appearing in 
the approximate operator, can be obtained from a Taylor Series 
Expansion in two variables. If a Taylor Series Expansion is 
centered on nodal point (m), values for all other nodal points 
in the grid may be expressed in terms of the nodal value at 
(m), its partial derivatives, and the grid spacing. Values for 
the nodal derivatives specified by the boundary conditions 
are also expressed in terms of the nodal point (m) by using 
the respective derivatives of the Taylor Series Expansion. 

Consequent evaluation of the required nodal values in 
terms of nodal point (m) and substitution into relationship 


(2.45) yields the following relationship for the error of 


approximation 
- n,+4 
2 i a Pin 
Bt a BE, Cs aac.) (2246) 
1=1 OX OY 
n,=2,3, ,2 
C.=6,7,-++5 


Es=CONSTANT 


Zi 





For the purpose of this work the above series was truncated 
after the first non-zero term in order to obtain a bound on 
the error for the respective approximate biharmonic operator. 
After truncation, the error for all of the approximate bihar- 
monic operators developed in this work was found to be propor- 
tional to at least the square of the mesh size (h) and a sixth 
partial derivative with respect to x and y. 

With the error proportional to at least the square of the 
mesh size (h) used for the finite difference approximation, 
resulting solutions should rapidly converge to the true solu- 
tion as the finite difference mesh is refined. The true 
solution should be obtained as the mesh size tends toward zero. 

To obtain an error proportional to at least the square 
of the mesh size spacing in all cases required the use of 
polynomial (2.2) for the development of all biharmonic molecules 


with boundary conditions included. 
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Pirie RESULTS 

iMiewmolLeecules Lor the various finite difference operators 
developed by the author in this work are presented in Tables 
mepethrough (8) . 

Table (1) represents the molecule for the general alyeaie= 
monic operator evaluated for the center nodal point. This 
molecule is applicable for all nodal points in a mesh which are 
at least two nodal points removed from all boundaries. 

Table (2) represents the molecules for nodal points ie 
3) and (2,3) to be used to evaluate the biharmonic operator 
when boundary conditions specify the first partial with 
respect to y along an edge. Such a case would be the solution 
of a flat plate bending problem with a built-in edge restraint. 

Table (3) represents the molecules for nodal points (1,3), 
(2,3), and (1,2) to be used to evaluate the biharmonic operator 
when boundary conditions specify the first partial along two 
mutual perpendicular edges. Such a case might be the solution 
for the deflection of a flat plate under uniform pressure with 
all edges built-in. 

Table (4) represents the molecules for nodal points (1,3) 
and (2,3) to be used for the biharmonic operator when boundary 
conditions specified by the harmonic operator are required. 
Such a boundary condition might be the moment for a flat plate 
with an edge simply supported. 

Table (5) represents the molecules for nodal points (1,3), 
ive wemancdei 2) to be used to numerically evaluate the bihar- 


monic operator when boundary conditions specify the second 


Zo 





Partial with respect to x along one edge of a boundary and 

the second partial with respect to y along another edge. Such 
a condition would exist in solving a plate bending problem 
when all edges are Simply supported. 

Table (6) represents the molecules for nodal points (1,3), 
(2,3), and (1,2) to be used to evaluate the biharmonic operator 
when boundary conditions specify the first partial with respect 
to y along one boundary and the harmonic operator along an 
adjacent boundary. Such a case would be for a flat plate 
with one edge built-in and an adjacent edge simply supported. 

Tables (7) and (8) represent the molecules for nodal 
Pemmeo (1,0) and (2,0) to be used to numerically evaluate 
the boundary conditions governed by equations (2.37) and 
258). 

Below each molecule a common factor is given by which the 
tabulated entries of the molecule must be multiplied. This 
common factor was generated by converting the original compu- 
tational values which were in terms of repeating decimal 
fractions to rational fractions. 

Each molecule is presented in a form Similar to the grid 
(Figures 1, 2, and 3) from which it was developed. Consequently 
each constant in a molecule applies to the respective nodal 
deflection or boundary condition value in the grid for which 
it was developed. 

An estimate of the error of approximation is given below 


each molecule. 
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IV. DISCUSSION 


A ieee Ehy sO PERI VvVED FINITE DIFFERENCE OPERATORS 
Confidence in the integrity of the finite difference 
operators developed in this work was established by comparing 
the solutions for the deflection of a flat plate under uniform 
pressure by the finite difference method with the exact solu- 

tion given by Timoshenko (Ref. 6) and Roark (Ref. 7). 

Each problem studied was solved using the molecules 
developed in this work and general finite difference techniques. 
Mesh sizes of L/4, L/6, L/8, L/10, L/12, where L is the width 
of the square plate studied, were used in successive solutions 
to the problem. The solution for the maximum deflection, 


6 of the plate was used aS a comparison. The value obtained 


max’ 
from successive refinement of mesh spacing for each problem was 
plotted versus L/N@, where N is the number of divisions on one 
Side of the plate. The resulting plots were then extrapolated 
to obtain the solution for a mesh having an infinite number of 
nodal points or the condition where the mesh spacing approaches 
zero. This extrapolated value is taken as the final solution 
for each problem and 1S compared with the exact solution. The 
above plots also indicate the convergence tendency of the 
solutions obtained with the molecules developed. 

The first problem solved was one of a Square plate under 
uniform pressure with all edges built-in. The problem was 
solved using the molecules developed in this work, (Present 


Method), and by means of the standard biharmonic molecule 


found in the literature and fictitious points, (Classical 
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Method). The results obtained from these methods are pre- 
Sented in Table (9). Figure (4) is the plot of successive 
solutions and illustrates the convergence trend for each 
method. 

A brief description is given in Appendix A of the solution 
process for the problems studied in this work. The solution 
process given is the one for the Present Method. 

The second problem solved was a Square plate under uniform 
pressure with all edges simply Supported. The results for the 
maximum deflection of the plate obtained from both methods of 
solution are presented in Table (10). The plot of the results 
iema1ven in Figure (5). 

The third problem solved wasS a Square plate under uniform 
pressure with two opposite edges built-in and the other two 
edges simply supported. A solution using the Classical Method, 
although possible, was not performed. The results for the 
maximum deflection at the center of the plate are given in 
Table (11) and the convergence plot by Figure (6). 

The last problem solved was a Square plate under uniform 
pressure with two opposite edges simply supported, one edge 
built-in and the remaining edge free. The Classical Method 
is not feasible for the solution of this problem since the 
relationship for fictitious points outside the free edge is 
unknown. The use of a L/4 mesh spacing was also determined to 
be too coarse fora valid solution of the problem. The reason 
being that equations in addition to the governing equilibrium 


equation must be written to account for the boundary conditions 
48 





MESH 
SIZE 


SPACING 


h 
L/4 
L/6 
L/8 
L/10 
7 1.2 


0 


MESH 
SIZE 


SPACING 


h 
L/4 
L/6 
L/8 
L/10 
ityAly 


0 


a et 


eee eee 


NO. OF PRESENT CLASSICAL EXACT 
MESH METHOD METHOD SOLN. 
DIVISIONS REF. 6 
N lo PL? /D mole PL? /D Om PL*/p 
4 1.0455 1.7995 
6 1.1293 1.5343 
8 inlets 1.4245 
10 1.2095 1.3696 
ii 1.2258 esses 
co 1.2628 1.2688 263" 


Table 9 - Maximum Deflection of 


in On All Edges Under 


a Square Plate Built- 
Mn eEOrtmMeEressure 


NO. OF PRESENT CLASSICAL Ae: 
MESH METHOD METHOD SOLN. 
DIVISIONS REF. 6 
N me PL? /p x1073 PL? /D x1073 PL? /D 
4 3.6555 4.0283 
6 B06 405 4.0483 
8 5.9265 4.0547 
10 geo 79 Ae oS 
a2 3.9979 4.0590 
co 4.0570 4.0630 4.0567 


Table 10 - Maximum Deflection For Square Plate Supported 
On All Edges Under Uniform Pressure 
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A PRESENT METHOD 
O CLASSICAL METHOD 
> EXACT VALUE 
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FIGURE 4 - Convergence Plot For Plate With 
All Edges Built-in 
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FIGURE 5 - Convergence Plot For Plate With 
All Edges Simply Supported 
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SPACING 
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Table 11 


MESH 
SIZE 
SPACING 


h 
L/6 
L/8 
L/10 
L/12 


0 


NO Or 
MESH 
DIVISIONS 


N 


4 


10 


Ps oN 
METHOD 


3 pit yp 


On 
1.6563 
1.7625 
ee 
1.8530 


1.9094 


xl0- 
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een? 
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- Maximum Deflection For Square Plate With 
Two Opposite Edges Built-in And The Other 
Two Simply Supported Under Uniform Pressure 
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Table 12 - Maximum.Deflection For Square Plate With 

One Hage Built-in, One Edge Simply Supported, 
One Edge Free And One Edge Simply Supported 
Under Uniform Pressure 
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FIGURE 6 - Convergence Plot For Plate Built-in 


On Two Opposite Edges and Simply 
Supported On The Other Two Edges 
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along the free edge, since these conditions were not incorpor- 
ated into the finite difference biharmonic operator. The 
coarse grid given by a L/4 mesh does not provide sufficient 
nodal points to write equations to establish complete equil- 
ibrium over the entire plate. The requirement of additional 
equations for this particular problem is made clearer in the 
solution process in Appendix A. The results for the maximum 
deflection of the plate are given in Table (12), and the 


convergence plot is illustrated by Figure (7). 


1D CONCLUSIONS 

Sets of mathematical molecules for the numerical approx- 
imation by finite differences of the biharmonic operator with 
various boundary conditions included were developed, as well 
as molecules for other differential operators. One consistent 
form of an interpolating polynomial was used to approximate 
the required two dimensional surface for the construction of 
all finite difference operators. 

Solutions for various plate bending problems obtained by 
using the developed finite difference operators are in 
excellent agreement with the exact solutions found in the 
literature. 

The convergence plots illustrate a straight convergence 
rate without oscillations for mesh spacings greater than L/6 
for all problems where the biharmonic operator and boundary 


conditions may be combined into one molecule. In Figures (4) 


5D 





and (5), solutions obtained by the present method are compared 
to the classical solutions with fictitious points. In one 
case the rate of convergence of the method is better than the 
classical method, in the other it is not. However the present 
method always gave a lower bound which is not always true for 
the classical method as shown in the two figures. 

The major advantage of the molecules developed is that 
they may be used to solve boundary values problems without the 
use of fictitious nodal points. Consequently, solutions to 
problems such as the plate with free edge boundary condition 
can now be solved uSing finite difference operators. 

One final conclusion concerning the overall aspect of 
the technique proposed should be made. Although the technique 
for developing finite difference operators proposed in this 
work has certain advantages over the classical method in the 
literature, it should be realised that other numerical methods 
exist for solving boundary value problems. The FINITE ELEMENT 
method has been applied very successfully and is covered in the 
literature very thoroughly, for example Refs. (8, 9 and 10). 

The method of finite differences using the consistent 
operators developed in this work is not superior to that of 
finite element methods for the solution of problems that may 
be formulated for a finite element sche:ie of solution, such as 
plate bending problems. However, difference operators may be 
developed, using the technique proposed in this work, to 


obtain an approximate solution to any problem for which the 
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governing equation is known. Thus they are applicable in the 
solution of problems that are not easily formulated into a 


mete element scheme of solution. 


Sy) 





APPENDIX A 
Pirie oe Ree ELON SOF SOLUTION PROCESS 

The solution of the plate bending problems by means of 
the finite difference operators developed in this work was 
performed by solving the governing equation Viw=p /D where 
w=w(x,y) is the displacement component perpendicular to the 
middle plane (x,y) of the plate. 

For the plate with all edges built-in, a mesh spacing 
of L/4 was used for the first solution. The applicable mesh 


is given below in Figure (8). 





FIGURE 8 - L/4 Finite Difference Mesh For A 
Square Plate With All Edges Built- 
in Or Simply Supported 
Due to symmetry and the fact that w=0 along the boundaries, 
there are only three unknown values of deflections as shown 
in Figure (8). For each unknown nodal deflection, the 
governing equation for plate equilibrium is written with the 
biharmonic operator being approximated by the appropriate 
difference operator developed in this work. The resulting 
set of simultaneous equations are then solved for the unknown 
deflections. For this problem the following molecules were 


used to express the governing equation 
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i. emcee wlenrrom Table (3) for co-ordinate (1,3) 
applied at nodal point one. 


ig Motecule from Table (3) for co-ordinate (2,3) 
applied at nodal point two. 

ce Molecule from Table (1) applied at nodal point 
three. 


The derivatives required for the nodal points along the bound- 
aries are set to zero, thus establishing zero slope along the 
built-in edges. 

For the plate with all edges Simply supported, the same 
mesh (Fig. (8)) is applicable for the first finite difference 
solution. 

The three unknown deflections are determined in this case 


by evaluating the governing equation with the following 


molecules 
ey Molecule from Table (5) for co-ordinate (1,3) 
applied at nodal point one. 
2s Molecule from Table (5) for co-ordinate (2,3) 
applied at nodal point two. 
Oy: Molecule from Table (1) applied at nodal point 
three. 


The derivatives required for the nodal points along the 
boundaries are set to zero, establishing zero moment for the 
Simply supported restraints. 

For a plate with two opposite edges built-in and two 
edges simply supported, the L/4 mesh shown in Fig. (9) is 


applicable 


ao 








FIGURE 9 - L/4 Finite Difference Mesh For A 
Square Plate With Top And Bottom 
Edges Built-in And Side Edges 
Simply Supported 
In this case, after applying symmetry and the fact of zero 
deflection along the edges, there are four unknowns to be 


determined. The required equations are written using the 


following molecules 


ale Molecule from Table (6) for co-ordinate (1,3) 
applied at nodal point one. 

Zs Molecule from Table (6) for co-ordinate (2,3) 
applied at nodal point two. 

Cie Molecule from Table (6) for co-ordinate (1,2) 
applied at nodal point three. 

4. Molecule from Table (1) applied at nodal point 
fours | 


As before, the values of the derivatives specified by the 
boundary conditions are set to zero. 

For the plate with built-in, simply supported, free and 
Simply supported edges, a mesh spacing of L/6 was used for 


the first finite difference solution as shown in Fig. (10). 
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FIGURE 10 - L/6 Finite Difference Mesh For A 

Square Plate With Top Edge Built- 

in, Sides Simply Supported, And 

Bottom Edge Free 
After applying symmetry and the fact that the deflection is 
zero on all boundaries except the free one, eighteen unknown 
deflections remain in the L/6 mesh. These unknowns are 
determined by solving the governing equation for nodal points 
1 through 12 and the boundary condition @€quations of shear 
pmaeiaoment: along a free edge for nodal points 16, 17, and 18. 
This requirement of uSing the nodal points along the free 
edge to specify the applicable boundary conditions at this 
location reduces the number of nodal points for which the 
governing equilibrium equation may be written. Consequently, 
an L/4 grid spacing was found to be too coarse for a valid 
solution of the problem. 


The required eighteen equations are written using the 


following finite difference molecules; 


ere Molecule from Table (6) for co-ordinate (1, 
3) applied at nodal point one. 

Dis Molecule from Table (6) for co-ordinate (2, 
3) applied at nodal point two. 

ry Molecule from Table (2) for co-ordinate (2, 


3) applied at nodal point three. 
Gal 





A Norecule=irom Table (6) £or co-ordinate (1, 
2) applied at nodal point four. 
5 OO 4 Molecule from Table (1) applied at nodal 


oral, pobmeses,0,6,9,L1, and 12. 
naz 
oO. Molecule from Table (4) for co-ordinate (2, 


3) applied at nodal points seven and ten. 
timetaromecase, the molecule ws laid on the 
finite difference mesh in the y direction 
so that the specified boundary condition 

constructed in the operator is along the 

Simply supported edge. 


Me. Molecule from Table (7) for co-ordinate (1, 
0) is applied at nodal point sixteen. 

AS ls < Molecule from Table (7) for co-ordinate (2, 
0) is applied at nodal points seventeen and 
eighteen. 

Gr Molecule from Table (8) for co-ordinate (1, 
0) is applied at nodal point sixteen. 

eye es. Molecule from Table (8) for co-ordinate (2, 
0) is applied at nodal points seventeen and 
eighteen. 


Values for the derivatives along the built-in and simply 
Supported edges required in the solution of this problem are 
set to zero to account for the specified boundary conditions. 
In order to obtain sufficient values for the deflection 
Beeomoarticular nodal point for the construction of a con- 
vergence plot, successive solutions are performed using a 


refinement in the mesh spacing for each solution. 
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